#setwd("C:/Users/erussek/forage_jsp/analysis")
setwd("/Users/evanrussek/forage_jsp/analysis")
rm(list = ls())
library(tidyr)
Attaching package: ‘tidyr’
The following object is masked from ‘package:Matrix’:
expand
library(dplyr)
Attaching package: ‘dplyr’
The following objects are masked from ‘package:stats’:
filter, lag
The following objects are masked from ‘package:base’:
intersect, setdiff, setequal, union
library(zoo)
Attaching package: ‘zoo’
The following objects are masked from ‘package:base’:
as.Date, as.Date.numeric
library(ggplot2)
library(ggpubr)
Loading required package: magrittr
Attaching package: ‘magrittr’
The following object is masked from ‘package:tidyr’:
extract
library(RcppRoll)
library(knitr)
data <- read.csv('data/run4_data.csv')
head(data)
data <- data %>% ungroup() %>% bind_cols(s_num = group_indices(., subjectID))
data <- data %>% mutate(subj = factor(s_num))
travel_keys = data %>% ungroup() %>% select(travel_key) %>% unique()
clean_subj_data <- function(subj_data, travel_keys){
travel_key_hard = travel_keys$travel_key[1]
travel_key_easy = travel_keys$travel_key[2]
subj_data <- data %>%select(round, phase, reward_obs, reward_true, lag, exit, start_reward, n_travel_steps,
travel_key, subjectID, trial_num, s_num, correct_key) %>% group_by(trial_num) %>%
mutate(press_num = 1:n(), round = as.factor(round),
rep_number = case_when(trial_num < 7 ~ 1, TRUE ~ 2)) %>%
mutate(phase = replace(phase, phase == "Harvest", "HARVEST"))
subj_data <- subj_data %>% group_by(trial_num) %>%
mutate(travel_key = replace(travel_key, travel_key == "", first(travel_key[travel_key != ""]))) %>%
mutate(travel_key_cond = case_when(travel_key == travel_key_hard ~ "HARD",
travel_key == travel_key_easy ~ "EASY")) %>%
filter(!is.na(travel_key_cond))
}
# clean all the data and bind
datalist <- list()
for (i in 1:10){
subj_data <- data %>% filter(s_num == i)
subj_data <- clean_subj_data(subj_data, travel_keys)
subj_data <- ungroup(subj_data)
datalist[[i]] <- subj_data
}
Unequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorUnequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorUnequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorUnequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorUnequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorUnequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorUnequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorUnequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorUnequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorUnequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vector
#cdata <- do.call(rbind,datalist)
cdata <- bind_rows(datalist)
# make a function to plot reward on each game...
plot_subj_reward_v_press <- function(subj_data){
s_num <- subj_data$s_num[1]
# reward plots which show the threshold...
rew_plot <- ggplot(subj_data, aes(x = press_num, y = reward_obs, group = round)) +
#geom_rect(data = t1_data, aes(xmin = press_num - 0.5, xmax = press_num + 0.5, ymin = -Inf, ymax = Inf, fill = round)) +
geom_point(aes(color = phase)) + facet_grid(n_travel_steps ~ travel_key_cond) + theme(legend.position = "none") + ylim(0,110) + ggtitle(paste("subj: ", s_num)) + ylab('reward collected') + xlab('button press number') + theme_minimal()
# plot_an <- annotate_figure(plot, top = paste('subj: ', s_num))
plot(rew_plot)
}
for (s in 1:22){
subj_data <- cdata %>% filter(s_num == s)
subj_data$round_num <- as.integer(as.character(subj_data$round))
#subj_data <- subj_data %>% group_by(travel_key_cond, n_travel_steps) %>% mutate(press_num = 1:n()) %>% ungroup()
plot_subj_reward_v_press(subj_data)
}






















subj_data$round_num
[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[41] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[81] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[121] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[161] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[201] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[241] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[281] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[321] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[361] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[401] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[441] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[481] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[521] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[561] 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[601] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[641] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[681] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[721] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[761] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
[801] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
[841] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
[881] 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4 4
[921] 4 4 4 4 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5
[961] 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 5 6 6 6
[ reached getOption("max.print") -- omitted 32240 entries ]
get_trial_exits <- function(thdata){
# get harvest
last_phase <- last(thdata$phase)
# go through data.. if last phase was harvest, remove that round...
if (last_phase == "HARVEST"){
last_round <- last(thdata$round)
# get the last reward ops...
last_reward_ob <- last(thdata$reward_obs[!is.na(thdata$reward_obs)])
#if (last_reward_ob > 8){
thdata <- thdata %>% filter(round != last_round)
#}
}
# now select harvest data
thdata <- thdata %>% filter(phase == "HARVEST", !is.na(reward_obs)) # find out why reward_true has so many .na
## get either last true reward observed
return_tbl <- thdata %>%
group_by(round) %>%
summarise(s_num = first(s_num), last_reward = last(reward_obs),
trial_num = first(trial_num),
n_travel_steps = first(n_travel_steps),
travel_key_cond = first(travel_key_cond),
rep_number = first(rep_number)) %>% ungroup()
return(return_tbl)
}
#sdata %>% group_by(trial_num) %>%
# do(get_trial_exits(.))
exit_data <- cdata %>% group_by(s_num, trial_num) %>%
do(get_trial_exits(.)) %>% ungroup() %>% mutate(subj = as.factor(s_num))
library(ggpubr)
exit_data <- exit_data %>% mutate(round_num = as.integer(as.character(round)))
for (s in 1:22){
p <- ggplot(exit_data %>% filter(s_num == s), aes(x = round_num, y = last_reward)) + geom_point() + geom_line() +
facet_grid(n_travel_steps ~ travel_key_cond) + theme_minimal() + ylim(0,100) + ylab('last reward collected') + xlab('harvest number')
#plot(p)
#ggarrange(p)
a <- annotate_figure(p, right = "# travel steps", top = "travel step difficulty", fig.lab = paste('subj: ' , s))
plot(a)
}
geom_path: Each group consists of only one observation. Do you need to adjust the group aesthetic?























# do some aggregating over exit thresholds...
# just get the mean for each subject for each timepoint, collapse accross rep number...
mn_exit <- exit_data %>% group_by(s_num, n_travel_steps, travel_key_cond) %>%
summarise(rep_exit_thresh = mean(last_reward), trial_num = mean(trial_num)) %>%
group_by(s_num, n_travel_steps, travel_key_cond) %>%
summarise(exit_thresh = mean(rep_exit_thresh), trial_num = mean(trial_num)) %>%
mutate(subj = as.factor(s_num)) %>% ungroup()
Unequal factor levels: coercing to characterbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vectorbinding character and factor vector, coercing into character vector
# now just plot this for each subject -- draw lines to show easy - hard effect
ggplot(mn_exit, aes(x = travel_key_cond, y = exit_thresh)) + geom_dotplot(binaxis = 'y', aes(fill = subj), dotsize = 1.2) + geom_line(aes(group = subj, color = subj), size = 1.2) +
facet_grid(.~n_travel_steps) + theme_minimal() + labs(y = 'Mean Last Reward Collected', x = 'Travel Sequence', subtitle = '# Travel Presses') #theme(legend.position = "none")
ggsave('reward_press2.png')
Saving 7.29 x 4.51 in image

ggplot(mn_exit, aes(x = factor(n_travel_steps), y = exit_thresh)) +
geom_dotplot(binaxis = 'y', aes(fill = subj), dotsize = 1.2) +
geom_line(aes(group = subj, color = subj), size = 1.2) +
facet_grid(.~travel_key_cond) +
theme_minimal() +
labs(y = 'Mean Last Reward Collected', x = 'Travel Key Condition', subtitle = '# Travel Presses')
#theme(legend.position = "none")
ggsave('reward_press.png')
Saving 7.29 x 4.51 in image

# plot the means for each of these...
gmn_exit <- mn_exit %>% group_by(n_travel_steps, travel_key_cond) %>% summarise(gm_thresh = mean(exit_thresh), gsd_thresh = sd(exit_thresh)/sqrt(n()))
ggplot(gmn_exit, aes(x = factor(n_travel_steps), y = gm_thresh, color = travel_key_cond)) +
geom_line(aes(group = travel_key_cond), size = 2) +
geom_errorbar(aes(ymin = gm_thresh - gsd_thresh, ymax = gm_thresh+gsd_thresh), width = .1, size = 2) +
geom_point(size = 4) + ylab('group mean exit threshold') + xlab('# travel steps') + labs(color = 'travel effort') + theme_minimal()
ggsave('group mean exit thresh.png')
Saving 7.29 x 4.51 in image

library(optimx)
library(lmerTest)
exit_model <- lmer(last_reward ~ n_travel_steps + travel_key_cond + (1 + n_travel_steps + travel_key_cond |subj), data = exit_data, control = lmerControl(optimizer ="optimx",optCtrl=list(method='nlminb')))
summary(exit_model)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: last_reward ~ n_travel_steps + travel_key_cond + (1 + n_travel_steps + travel_key_cond | subj)
Data: exit_data
Control: lmerControl(optimizer = "optimx", optCtrl = list(method = "nlminb"))
REML criterion at convergence: 6956.1
Scaled residuals:
Min 1Q Median 3Q Max
-4.1827 -0.5108 0.0002 0.5046 5.3563
Random effects:
Groups Name Variance Std.Dev. Corr
subj (Intercept) 290.72330 17.0506
n_travel_steps 0.06055 0.2461 -0.27
travel_key_condHARD 68.01256 8.2470 -0.25 0.26
Residual 87.82008 9.3712
Number of obs: 928, groups: subj, 22
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 54.95291 3.71001 20.77742 14.812 1.63e-12 ***
n_travel_steps -0.24696 0.06318 19.11207 -3.909 0.000935 ***
travel_key_condHARD -5.06293 1.89212 17.91918 -2.676 0.015465 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) n_trv_
n_trvl_stps -0.309
trvl_k_HARD -0.256 0.206
library(broom.mixed)
td <- tidy(exit_model, conf.int = TRUE)
tdsome <- td %>% filter(term == 'n_travel_steps' | term == 'travel_key_condHARD')
tdsome <- tdsome %>% mutate(parameter = case_when(term == 'n_travel_steps' ~ 'travel step', term == 'travel_key_condHARD' ~ 'travel effort cond'))
ggplot(tdsome, aes(estimate,parameter,color = parameter)) +
geom_point(size = 2) +
geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), size = 1.2) + theme(legend.position = "none") + theme_minimal() + ylab('fixed effect') + xlab('estimate (# points)') + ggtitle('linear mixed effects model - effect on exit threshold') + theme(axis.text.y = element_text(angle = 90, hjust = .5))
ggsave('exit model.png')
Saving 7.29 x 4.51 in image

# reaction times...
## a bit more sensible...
lag_data <- cdata %>% select(s_num, rep_number, travel_key_cond, n_travel_steps, phase, trial_num, lag, correct_key, round, trial_num) %>% mutate(log_lag = log(lag)) %>% group_by(s_num, trial_num, round, phase) %>% slice(-1) %>% ungroup() %>% mutate(subj = as.factor(s_num))
ggplot(lag_data, aes(x=lag)) + geom_histogram(binwidth = 10) + xlim(50,300) + facet_wrap( ~ s_num)

#ggplot(lag_data, aes(x=log_lag)) + geom_histogram(binwidth = .01) + facet_wrap( ~ s_num)
filt_lag <- lag_data %>% group_by(s_num) %>% filter(log_lag < (mean(log_lag) + 1.5*sd(log_lag)) , log_lag > (mean(log_lag) - 1.5*sd(log_lag))) %>% ungroup()
ggplot(filt_lag, aes(x=lag)) + geom_histogram(binwidth = 10) + facet_wrap(~s_num) + xlim(50,350)

ggplot(filt_lag, aes(x=log_lag)) + geom_histogram(binwidth = .02) + facet_wrap(~s_num)

## with a run, plot the mean lag
round_filt_lag <- filt_lag %>%
group_by(s_num, n_travel_steps, travel_key_cond, round, phase) %>%
summarise(trial_num = first(trial_num), mean_lag = median(lag), mean_log_lag = median(log_lag)) %>% mutate(round_num = as.integer(as.character(round)))
for (s in 1:22){
# plot mean lag as a function of round number
p <- ggplot(round_filt_lag %>% filter(s_num == s), aes(x = round_num, y = mean_log_lag, color = phase)) + geom_point(size = 2) +facet_grid(n_travel_steps ~ travel_key_cond) + geom_line(size = 1.2) + xlab('harvest/travel number') + ylab('mean log lag') + ggtitle(paste('subj: ', s)) + theme_minimal()
plot(p)
}






















NA
NA
library(plotrix)
cond_round_filt_lag <- round_filt_lag %>%
group_by(n_travel_steps, travel_key_cond, round_num, phase) %>%
summarise(mean_log_lag = median(mean_log_lag), sd_lag = std.error(mean_log_lag))
p <- ggplot(cond_round_filt_lag, aes(x = round_num, y = mean_log_lag, color = phase)) + geom_point(size = 2) +facet_grid(n_travel_steps ~ travel_key_cond) + geom_line(size = 1.2) + geom_linerange(aes(ymin = mean_log_lag - sd_lag, ymax = mean_log_lag + sd_lag)) +
xlab('harvest/travel number') + ylab('mean log lag') + ggtitle('group: ') + theme_minimal()
plot(p)

cond_round_filt_lag
trial_filt_lag <- round_filt_lag %>%
group_by(s_num, n_travel_steps, travel_key_cond, phase) %>%
summarise(trial_num = first(trial_num),
mean_lag = mean(mean_lag),
mean_log_lag = mean(mean_log_lag)) %>% ungroup() %>% mutate(s_num = as.factor(s_num))
ggplot(trial_filt_lag,
aes(x = factor(n_travel_steps), y = mean_log_lag, fill = s_num)) +
geom_point(aes(color = s_num), size = 2)+
geom_line(aes(group = s_num, color = s_num), size = 1.2) +
facet_grid(phase ~ travel_key_cond) + theme_minimal() + xlab('# travel steps') + ylab('log lag time')

ggplot(trial_filt_lag,
aes(x = travel_key_cond, y = mean_lag, fill = s_num)) +
geom_point(aes(color = s_num), size = 2)+
geom_line(aes(group = s_num, color = s_num), size = 1.2) +
facet_grid(phase ~ n_travel_steps) + theme_minimal()

group_lag <- trial_filt_lag %>% group_by(n_travel_steps, travel_key_cond, phase) %>% summarise(gm_lag = mean(mean_lag), gml_lag = mean(mean_log_lag)) %>% ungroup()
ggplot(group_lag, aes(x = factor(n_travel_steps), y = gml_lag, color = travel_key_cond)) +
geom_line(aes(group = travel_key_cond), size = 2) + facet_wrap(~phase) + geom_point(size = 2) +xlab('# travel steps') +
ylab('group mean log lag') + theme_minimal()

NA
NA
NA
harvest_lag <- round_filt_lag %>% ungroup() %>% filter(phase == "HARVEST") %>% mutate(subj = as.factor(s_num))
head(harvest_lag)
hl_model <- lmer(scale(mean_log_lag) ~ n_travel_steps + travel_key_cond + round_num + (1 + n_travel_steps+ travel_key_cond + round_num |subj), data = harvest_lag, control = lmerControl(optimizer ="optimx",optCtrl=list(method='nlminb')))
summary(hl_model)
Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: scale(mean_log_lag) ~ n_travel_steps + travel_key_cond + round_num +
(1 + n_travel_steps + travel_key_cond + round_num | subj)
Data: harvest_lag
Control: lmerControl(optimizer = "optimx", optCtrl = list(method = "nlminb"))
REML criterion at convergence: 1236.7
Scaled residuals:
Min 1Q Median 3Q Max
-3.6966 -0.4726 -0.0215 0.4743 7.2375
Random effects:
Groups Name Variance Std.Dev. Corr
subj (Intercept) 0.7512496 0.86675
n_travel_steps 0.0001744 0.01320 -0.06
travel_key_condHARD 0.1057507 0.32519 -0.08 -0.28
round_num 0.0001980 0.01407 -0.33 0.27 -0.07
Residual 0.1563930 0.39547
Number of obs: 1029, groups: subj, 22
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -0.067679 0.189199 21.299098 -0.358 0.724080
n_travel_steps -0.000230 0.003159 21.798357 -0.073 0.942619
travel_key_condHARD 0.076090 0.074627 19.314980 1.020 0.320517
round_num 0.023682 0.005448 15.361270 4.347 0.000546 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Correlation of Fixed Effects:
(Intr) n_trv_ t__HAR
n_trvl_stps -0.128
trvl_k_HARD -0.103 -0.220
round_num -0.298 0.228 0.012
library(broom.mixed)
td <- tidy(hl_model, conf.int = TRUE)
tdsome <- td %>% filter(term == 'n_travel_steps' | term == 'travel_key_condHARD' | term == 'round_num')
tdsome <- tdsome %>% mutate(parameter = case_when(term == 'n_travel_steps' ~ 'travel step', term == 'travel_key_condHARD' ~ 'travel effort cond', term == 'round_num' ~ 'round number'))
ggplot(tdsome, aes(estimate,parameter,color = parameter)) +
geom_point(size = 2) +
geom_errorbarh(aes(xmin = conf.low, xmax = conf.high), size = 1.2) + theme(legend.position = "none") + theme_minimal() + ylab('fixed effect') + xlab('estimate (# points)') + ggtitle('linear mixed effects model - effect on lag') + theme(axis.text.y = element_text(angle = 90, hjust = .5))

LS0tCnRpdGxlOiAiUiBOb3RlYm9vayIKb3V0cHV0OgogIGh0bWxfbm90ZWJvb2s6IGRlZmF1bHQKICBwZGZfZG9jdW1lbnQ6IGRlZmF1bHQKLS0tCgoKCmBgYHtyfQojc2V0d2QoIkM6L1VzZXJzL2VydXNzZWsvZm9yYWdlX2pzcC9hbmFseXNpcyIpCnNldHdkKCIvVXNlcnMvZXZhbnJ1c3Nlay9mb3JhZ2VfanNwL2FuYWx5c2lzIikKcm0obGlzdCA9IGxzKCkpCmxpYnJhcnkodGlkeXIpCmxpYnJhcnkoZHBseXIpCmxpYnJhcnkoem9vKQpsaWJyYXJ5KGdncGxvdDIpCmxpYnJhcnkoZ2dwdWJyKQpsaWJyYXJ5KFJjcHBSb2xsKQpsaWJyYXJ5KGtuaXRyKQpgYGAKCmBgYHtyfQpkYXRhIDwtIHJlYWQuY3N2KCdkYXRhL3J1bjRfZGF0YS5jc3YnKQpoZWFkKGRhdGEpCmRhdGEgPC0gZGF0YSAlPiUgdW5ncm91cCgpICU+JSBiaW5kX2NvbHMoc19udW0gPSBncm91cF9pbmRpY2VzKC4sIHN1YmplY3RJRCkpCmRhdGEgPC0gZGF0YSAlPiUgbXV0YXRlKHN1YmogPSBmYWN0b3Ioc19udW0pKQpgYGAKCmBgYHtyfQp0cmF2ZWxfa2V5cyA9IGRhdGEgJT4lIHVuZ3JvdXAoKSAlPiUgc2VsZWN0KHRyYXZlbF9rZXkpICU+JSB1bmlxdWUoKQpjbGVhbl9zdWJqX2RhdGEgPC0gZnVuY3Rpb24oc3Vial9kYXRhLCB0cmF2ZWxfa2V5cyl7CiAgdHJhdmVsX2tleV9oYXJkID0gdHJhdmVsX2tleXMkdHJhdmVsX2tleVsxXQogIHRyYXZlbF9rZXlfZWFzeSA9IHRyYXZlbF9rZXlzJHRyYXZlbF9rZXlbMl0KICAKICBzdWJqX2RhdGEgPC0gZGF0YSAlPiVzZWxlY3Qocm91bmQsIHBoYXNlLCByZXdhcmRfb2JzLCByZXdhcmRfdHJ1ZSwgbGFnLCBleGl0LCBzdGFydF9yZXdhcmQsIG5fdHJhdmVsX3N0ZXBzLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB0cmF2ZWxfa2V5LCBzdWJqZWN0SUQsIHRyaWFsX251bSwgc19udW0sIGNvcnJlY3Rfa2V5KSAlPiUgZ3JvdXBfYnkodHJpYWxfbnVtKSAlPiUKICAgIG11dGF0ZShwcmVzc19udW0gPSAxOm4oKSwgcm91bmQgPSBhcy5mYWN0b3Iocm91bmQpLCAKICAgICAgICAgICByZXBfbnVtYmVyID0gY2FzZV93aGVuKHRyaWFsX251bSA8IDcgfiAxLCBUUlVFIH4gMikpICU+JQogICAgbXV0YXRlKHBoYXNlID0gcmVwbGFjZShwaGFzZSwgcGhhc2UgPT0gIkhhcnZlc3QiLCAiSEFSVkVTVCIpKQoKICBzdWJqX2RhdGEgPC0gc3Vial9kYXRhICU+JSBncm91cF9ieSh0cmlhbF9udW0pICU+JSAKICAgIG11dGF0ZSh0cmF2ZWxfa2V5ID0gcmVwbGFjZSh0cmF2ZWxfa2V5LCB0cmF2ZWxfa2V5ID09ICIiLCBmaXJzdCh0cmF2ZWxfa2V5W3RyYXZlbF9rZXkgIT0gIiJdKSkpICU+JQogICAgbXV0YXRlKHRyYXZlbF9rZXlfY29uZCA9IGNhc2Vfd2hlbih0cmF2ZWxfa2V5ID09IHRyYXZlbF9rZXlfaGFyZCAgIH4gIkhBUkQiLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB0cmF2ZWxfa2V5ID09IHRyYXZlbF9rZXlfZWFzeSAgfiAiRUFTWSIpKSAlPiUKICAgIGZpbHRlcighaXMubmEodHJhdmVsX2tleV9jb25kKSkKfQoKIyBjbGVhbiBhbGwgdGhlIGRhdGEgYW5kIGJpbmQKZGF0YWxpc3QgPC0gbGlzdCgpCmZvciAoaSBpbiAxOjEwKXsKICBzdWJqX2RhdGEgPC0gZGF0YSAlPiUgZmlsdGVyKHNfbnVtID09IGkpCiAgc3Vial9kYXRhIDwtIGNsZWFuX3N1YmpfZGF0YShzdWJqX2RhdGEsIHRyYXZlbF9rZXlzKQogIHN1YmpfZGF0YSA8LSB1bmdyb3VwKHN1YmpfZGF0YSkKICBkYXRhbGlzdFtbaV1dIDwtIHN1YmpfZGF0YQp9CiNjZGF0YSA8LSBkby5jYWxsKHJiaW5kLGRhdGFsaXN0KQpjZGF0YSA8LSBiaW5kX3Jvd3MoZGF0YWxpc3QpCmBgYApgYGB7cn0KIyBtYWtlIGEgZnVuY3Rpb24gdG8gcGxvdCByZXdhcmQgb24gZWFjaCBnYW1lLi4uIApwbG90X3N1YmpfcmV3YXJkX3ZfcHJlc3MgPC0gZnVuY3Rpb24oc3Vial9kYXRhKXsKICAKICBzX251bSA8LSBzdWJqX2RhdGEkc19udW1bMV0KICAKICAjIHJld2FyZCBwbG90cyB3aGljaCBzaG93IHRoZSB0aHJlc2hvbGQuLi4gCiAgcmV3X3Bsb3QgPC0gZ2dwbG90KHN1YmpfZGF0YSwgYWVzKHggPSBwcmVzc19udW0sIHkgPSByZXdhcmRfb2JzLCBncm91cCA9IHJvdW5kKSkgICsKICAgICNnZW9tX3JlY3QoZGF0YSA9IHQxX2RhdGEsIGFlcyh4bWluID0gcHJlc3NfbnVtIC0gMC41LCB4bWF4ID0gcHJlc3NfbnVtICsgMC41LCB5bWluID0gLUluZiwgeW1heCA9IEluZiwgZmlsbCA9IHJvdW5kKSkgKwogICAgZ2VvbV9wb2ludChhZXMoY29sb3IgPSBwaGFzZSkpICsgZmFjZXRfZ3JpZChuX3RyYXZlbF9zdGVwcyB+IHRyYXZlbF9rZXlfY29uZCkgKyB0aGVtZShsZWdlbmQucG9zaXRpb24gPSAibm9uZSIpICsgeWxpbSgwLDExMCkgKyBnZ3RpdGxlKHBhc3RlKCJzdWJqOiAiLCBzX251bSkpICsgeWxhYigncmV3YXJkIGNvbGxlY3RlZCcpICsgeGxhYignYnV0dG9uIHByZXNzIG51bWJlcicpICsgdGhlbWVfbWluaW1hbCgpCiAgCiAjIHBsb3RfYW4gPC0gYW5ub3RhdGVfZmlndXJlKHBsb3QsIHRvcCA9IHBhc3RlKCdzdWJqOiAnLCBzX251bSkpCiAgcGxvdChyZXdfcGxvdCkKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAKfQpmb3IgKHMgaW4gMToyMil7CiAgc3Vial9kYXRhIDwtIGNkYXRhICU+JSBmaWx0ZXIoc19udW0gPT0gcykKICBzdWJqX2RhdGEkcm91bmRfbnVtIDwtIGFzLmludGVnZXIoYXMuY2hhcmFjdGVyKHN1YmpfZGF0YSRyb3VuZCkpCgogICNzdWJqX2RhdGEgPC0gc3Vial9kYXRhICU+JSBncm91cF9ieSh0cmF2ZWxfa2V5X2NvbmQsIG5fdHJhdmVsX3N0ZXBzKSAlPiUgbXV0YXRlKHByZXNzX251bSA9IDE6bigpKSAlPiUgdW5ncm91cCgpCiAgcGxvdF9zdWJqX3Jld2FyZF92X3ByZXNzKHN1YmpfZGF0YSkgIAp9CmBgYApgYGB7cn0Kc3Vial9kYXRhJHJvdW5kX251bSAKYGBgCgoKYGBge3J9CmdldF90cmlhbF9leGl0cyA8LSBmdW5jdGlvbih0aGRhdGEpewogICAgIyBnZXQgaGFydmVzdAogIGxhc3RfcGhhc2UgPC0gbGFzdCh0aGRhdGEkcGhhc2UpCiAgCiAgIyBnbyB0aHJvdWdoIGRhdGEuLiBpZiBsYXN0IHBoYXNlIHdhcyBoYXJ2ZXN0LCByZW1vdmUgdGhhdCByb3VuZC4uLiAKICBpZiAobGFzdF9waGFzZSA9PSAiSEFSVkVTVCIpewogICAgbGFzdF9yb3VuZCA8LSBsYXN0KHRoZGF0YSRyb3VuZCkKICAgICMgZ2V0IHRoZSBsYXN0IHJld2FyZCBvcHMuLi4KICAgIGxhc3RfcmV3YXJkX29iIDwtIGxhc3QodGhkYXRhJHJld2FyZF9vYnNbIWlzLm5hKHRoZGF0YSRyZXdhcmRfb2JzKV0pCiAgICAKICAgICNpZiAobGFzdF9yZXdhcmRfb2IgPiA4KXsKICAgIHRoZGF0YSA8LSB0aGRhdGEgJT4lIGZpbHRlcihyb3VuZCAhPSBsYXN0X3JvdW5kKQogICAgI30KICB9CiAgCiAgIyBub3cgc2VsZWN0IGhhcnZlc3QgZGF0YQogIHRoZGF0YSA8LSB0aGRhdGEgJT4lIGZpbHRlcihwaGFzZSA9PSAiSEFSVkVTVCIsICFpcy5uYShyZXdhcmRfb2JzKSkgIyBmaW5kIG91dCB3aHkgcmV3YXJkX3RydWUgaGFzIHNvIG1hbnkgLm5hCiAgCiAgIyMgZ2V0IGVpdGhlciBsYXN0IHRydWUgcmV3YXJkIG9ic2VydmVkCiAgcmV0dXJuX3RibCA8LSB0aGRhdGEgJT4lIAogICAgICAgICAgICAgICAgICBncm91cF9ieShyb3VuZCkgJT4lIAogICAgICAgICAgICAgICAgICAgICAgc3VtbWFyaXNlKHNfbnVtID0gZmlyc3Qoc19udW0pLCBsYXN0X3Jld2FyZCA9IGxhc3QocmV3YXJkX29icyksCiAgICAgICAgICAgICAgICAgICAgICAgICAgdHJpYWxfbnVtID0gZmlyc3QodHJpYWxfbnVtKSwKICAgICAgICAgICAgICAgICAgICAgICAgICBuX3RyYXZlbF9zdGVwcyA9IGZpcnN0KG5fdHJhdmVsX3N0ZXBzKSwKICAgICAgICAgICAgICAgICAgICAgICAgICB0cmF2ZWxfa2V5X2NvbmQgPSBmaXJzdCh0cmF2ZWxfa2V5X2NvbmQpLAogICAgICAgICAgICAgICAgICAgICAgICAgIHJlcF9udW1iZXIgPSBmaXJzdChyZXBfbnVtYmVyKSkgJT4lIHVuZ3JvdXAoKQogIHJldHVybihyZXR1cm5fdGJsKQp9Cgojc2RhdGEgJT4lIGdyb3VwX2J5KHRyaWFsX251bSkgJT4lCiMgIGRvKGdldF90cmlhbF9leGl0cyguKSkKCmV4aXRfZGF0YSA8LSBjZGF0YSAlPiUgZ3JvdXBfYnkoc19udW0sIHRyaWFsX251bSkgJT4lCiAgZG8oZ2V0X3RyaWFsX2V4aXRzKC4pKSAlPiUgdW5ncm91cCgpICU+JSBtdXRhdGUoc3ViaiA9IGFzLmZhY3RvcihzX251bSkpCmBgYAoKYGBge3J9CmxpYnJhcnkoZ2dwdWJyKQpleGl0X2RhdGEgPC0gZXhpdF9kYXRhICU+JSBtdXRhdGUocm91bmRfbnVtID0gYXMuaW50ZWdlcihhcy5jaGFyYWN0ZXIocm91bmQpKSkKZm9yIChzIGluIDE6MjIpewogIHAgPC0gZ2dwbG90KGV4aXRfZGF0YSAlPiUgZmlsdGVyKHNfbnVtID09IHMpLCBhZXMoeCA9IHJvdW5kX251bSwgeSA9IGxhc3RfcmV3YXJkKSkgKyBnZW9tX3BvaW50KCkgKyBnZW9tX2xpbmUoKSArCiAgZmFjZXRfZ3JpZChuX3RyYXZlbF9zdGVwcyB+IHRyYXZlbF9rZXlfY29uZCkgKyB0aGVtZV9taW5pbWFsKCkgICsgeWxpbSgwLDEwMCkgKyB5bGFiKCdsYXN0IHJld2FyZCBjb2xsZWN0ZWQnKSArIHhsYWIoJ2hhcnZlc3QgbnVtYmVyJykKICAKICAjcGxvdChwKQogICNnZ2FycmFuZ2UocCkKICBhIDwtIGFubm90YXRlX2ZpZ3VyZShwLCByaWdodCA9ICIjIHRyYXZlbCBzdGVwcyIsIHRvcCA9ICJ0cmF2ZWwgc3RlcCBkaWZmaWN1bHR5IiwgZmlnLmxhYiA9IHBhc3RlKCdzdWJqOiAnICwgcykpCiAgcGxvdChhKQp9CgpgYGAKCmBgYHtyfQojIGRvIHNvbWUgYWdncmVnYXRpbmcgb3ZlciBleGl0IHRocmVzaG9sZHMuLi4gCiMganVzdCBnZXQgdGhlIG1lYW4gZm9yIGVhY2ggc3ViamVjdCBmb3IgZWFjaCB0aW1lcG9pbnQsIGNvbGxhcHNlIGFjY3Jvc3MgcmVwIG51bWJlci4uLgptbl9leGl0IDwtIGV4aXRfZGF0YSAlPiUgZ3JvdXBfYnkoc19udW0sIG5fdHJhdmVsX3N0ZXBzLCB0cmF2ZWxfa2V5X2NvbmQpICU+JQogIHN1bW1hcmlzZShyZXBfZXhpdF90aHJlc2ggPSBtZWFuKGxhc3RfcmV3YXJkKSwgdHJpYWxfbnVtID0gbWVhbih0cmlhbF9udW0pKSAlPiUgCiAgZ3JvdXBfYnkoc19udW0sIG5fdHJhdmVsX3N0ZXBzLCB0cmF2ZWxfa2V5X2NvbmQpICU+JQogIHN1bW1hcmlzZShleGl0X3RocmVzaCA9IG1lYW4ocmVwX2V4aXRfdGhyZXNoKSwgdHJpYWxfbnVtID0gbWVhbih0cmlhbF9udW0pKSAlPiUKICBtdXRhdGUoc3ViaiA9IGFzLmZhY3RvcihzX251bSkpICU+JSB1bmdyb3VwKCkKCiMgbm93IGp1c3QgcGxvdCB0aGlzIGZvciBlYWNoIHN1YmplY3QgLS0gZHJhdyBsaW5lcyB0byBzaG93IGVhc3kgLSBoYXJkIGVmZmVjdApnZ3Bsb3QobW5fZXhpdCwgYWVzKHggPSB0cmF2ZWxfa2V5X2NvbmQsIHkgPSBleGl0X3RocmVzaCkpICsgZ2VvbV9kb3RwbG90KGJpbmF4aXMgPSAneScsIGFlcyhmaWxsID0gc3ViaiksIGRvdHNpemUgPSAxLjIpICsgZ2VvbV9saW5lKGFlcyhncm91cCA9IHN1YmosIGNvbG9yID0gc3ViaiksIHNpemUgPSAxLjIpICsKICBmYWNldF9ncmlkKC5+bl90cmF2ZWxfc3RlcHMpICsgdGhlbWVfbWluaW1hbCgpICsgbGFicyh5ID0gJ01lYW4gTGFzdCBSZXdhcmQgQ29sbGVjdGVkJywgeCA9ICdUcmF2ZWwgU2VxdWVuY2UnLCBzdWJ0aXRsZSA9ICcjIFRyYXZlbCBQcmVzc2VzJykgI3RoZW1lKGxlZ2VuZC5wb3NpdGlvbiA9ICJub25lIikKZ2dzYXZlKCdyZXdhcmRfcHJlc3MyLnBuZycpCgoKZ2dwbG90KG1uX2V4aXQsIGFlcyh4ID0gZmFjdG9yKG5fdHJhdmVsX3N0ZXBzKSwgeSA9IGV4aXRfdGhyZXNoKSkgKwogIGdlb21fZG90cGxvdChiaW5heGlzID0gJ3knLCBhZXMoZmlsbCA9IHN1YmopLCBkb3RzaXplID0gMS4yKSArIAogIGdlb21fbGluZShhZXMoZ3JvdXAgPSBzdWJqLCBjb2xvciA9IHN1YmopLCBzaXplID0gMS4yKSArCiAgZmFjZXRfZ3JpZCgufnRyYXZlbF9rZXlfY29uZCkgKyAKICB0aGVtZV9taW5pbWFsKCkgKyAKICBsYWJzKHkgPSAnTWVhbiBMYXN0IFJld2FyZCBDb2xsZWN0ZWQnLCB4ID0gJ1RyYXZlbCBLZXkgQ29uZGl0aW9uJywgc3VidGl0bGUgPSAnIyBUcmF2ZWwgUHJlc3NlcycpIAogICN0aGVtZShsZWdlbmQucG9zaXRpb24gPSAibm9uZSIpCmdnc2F2ZSgncmV3YXJkX3ByZXNzLnBuZycpCgpgYGAKYGBge3J9CiMgcGxvdCB0aGUgbWVhbnMgZm9yIGVhY2ggb2YgdGhlc2UuLi4KZ21uX2V4aXQgPC0gbW5fZXhpdCAlPiUgZ3JvdXBfYnkobl90cmF2ZWxfc3RlcHMsIHRyYXZlbF9rZXlfY29uZCkgJT4lIHN1bW1hcmlzZShnbV90aHJlc2ggPSBtZWFuKGV4aXRfdGhyZXNoKSwgZ3NkX3RocmVzaCA9IHNkKGV4aXRfdGhyZXNoKS9zcXJ0KG4oKSkpCgpnZ3Bsb3QoZ21uX2V4aXQsIGFlcyh4ID0gZmFjdG9yKG5fdHJhdmVsX3N0ZXBzKSwgeSA9IGdtX3RocmVzaCwgY29sb3IgPSB0cmF2ZWxfa2V5X2NvbmQpKSArCiAgZ2VvbV9saW5lKGFlcyhncm91cCA9IHRyYXZlbF9rZXlfY29uZCksIHNpemUgPSAyKSArCiAgZ2VvbV9lcnJvcmJhcihhZXMoeW1pbiA9IGdtX3RocmVzaCAtIGdzZF90aHJlc2gsIHltYXggPSBnbV90aHJlc2grZ3NkX3RocmVzaCksIHdpZHRoID0gLjEsIHNpemUgPSAyKSArIApnZW9tX3BvaW50KHNpemUgPSA0KSArIHlsYWIoJ2dyb3VwIG1lYW4gZXhpdCB0aHJlc2hvbGQnKSArIHhsYWIoJyMgdHJhdmVsIHN0ZXBzJykgKyBsYWJzKGNvbG9yID0gJ3RyYXZlbCBlZmZvcnQnKSArIHRoZW1lX21pbmltYWwoKQpnZ3NhdmUoJ2dyb3VwIG1lYW4gZXhpdCB0aHJlc2gucG5nJykKCmBgYAoKYGBge3J9CgpsaWJyYXJ5KG9wdGlteCkKbGlicmFyeShsbWVyVGVzdCkKZXhpdF9tb2RlbCA8LSBsbWVyKGxhc3RfcmV3YXJkIH4gbl90cmF2ZWxfc3RlcHMgKyB0cmF2ZWxfa2V5X2NvbmQgICsgKDEgICsgbl90cmF2ZWxfc3RlcHMgKyB0cmF2ZWxfa2V5X2NvbmQgfHN1YmopLCBkYXRhID0gZXhpdF9kYXRhLCAgY29udHJvbCA9IGxtZXJDb250cm9sKG9wdGltaXplciA9Im9wdGlteCIsb3B0Q3RybD1saXN0KG1ldGhvZD0nbmxtaW5iJykpKQpzdW1tYXJ5KGV4aXRfbW9kZWwpCiAgICAgICAgICAgICAgICAgICAKYGBgCgpgYGB7cn0KbGlicmFyeShicm9vbS5taXhlZCkKdGQgPC0gdGlkeShleGl0X21vZGVsLCBjb25mLmludCA9IFRSVUUpCnRkc29tZSA8LSB0ZCAlPiUgZmlsdGVyKHRlcm0gPT0gJ25fdHJhdmVsX3N0ZXBzJyB8IHRlcm0gPT0gJ3RyYXZlbF9rZXlfY29uZEhBUkQnKQoKdGRzb21lIDwtIHRkc29tZSAlPiUgbXV0YXRlKHBhcmFtZXRlciA9IGNhc2Vfd2hlbih0ZXJtID09ICduX3RyYXZlbF9zdGVwcycgfiAndHJhdmVsIHN0ZXAnLCB0ZXJtID09ICd0cmF2ZWxfa2V5X2NvbmRIQVJEJyB+ICd0cmF2ZWwgZWZmb3J0IGNvbmQnKSkKCmdncGxvdCh0ZHNvbWUsIGFlcyhlc3RpbWF0ZSxwYXJhbWV0ZXIsY29sb3IgPSBwYXJhbWV0ZXIpKSArIAogIGdlb21fcG9pbnQoc2l6ZSA9IDIpICArCmdlb21fZXJyb3JiYXJoKGFlcyh4bWluID0gY29uZi5sb3csIHhtYXggPSBjb25mLmhpZ2gpLCBzaXplID0gMS4yKSArIHRoZW1lKGxlZ2VuZC5wb3NpdGlvbiA9ICJub25lIikgKyB0aGVtZV9taW5pbWFsKCkgKyB5bGFiKCdmaXhlZCBlZmZlY3QnKSArIHhsYWIoJ2VzdGltYXRlICgjIHBvaW50cyknKSArIGdndGl0bGUoJ2xpbmVhciBtaXhlZCBlZmZlY3RzIG1vZGVsIC0gZWZmZWN0IG9uIGV4aXQgdGhyZXNob2xkJykgKyB0aGVtZShheGlzLnRleHQueSA9IGVsZW1lbnRfdGV4dChhbmdsZSA9IDkwLCBoanVzdCA9IC41KSkKZ2dzYXZlKCdleGl0IG1vZGVsLnBuZycpCmBgYAoKYGBge3J9CiMgcmVhY3Rpb24gdGltZXMuLi4gCgoKIyMgYSBiaXQgbW9yZSBzZW5zaWJsZS4uLiAKbGFnX2RhdGEgPC0gY2RhdGEgJT4lIHNlbGVjdChzX251bSwgcmVwX251bWJlciwgdHJhdmVsX2tleV9jb25kLCBuX3RyYXZlbF9zdGVwcywgcGhhc2UsIHRyaWFsX251bSwgbGFnLCBjb3JyZWN0X2tleSwgcm91bmQsIHRyaWFsX251bSkgJT4lIG11dGF0ZShsb2dfbGFnID0gbG9nKGxhZykpICU+JSBncm91cF9ieShzX251bSwgdHJpYWxfbnVtLCByb3VuZCwgcGhhc2UpICU+JSBzbGljZSgtMSkgJT4lIHVuZ3JvdXAoKSAlPiUgbXV0YXRlKHN1YmogPSBhcy5mYWN0b3Ioc19udW0pKQoKZ2dwbG90KGxhZ19kYXRhLCBhZXMoeD1sYWcpKSArIGdlb21faGlzdG9ncmFtKGJpbndpZHRoID0gMTApICsgeGxpbSg1MCwzMDApICsgZmFjZXRfd3JhcCggfiBzX251bSkKCiNnZ3Bsb3QobGFnX2RhdGEsIGFlcyh4PWxvZ19sYWcpKSAgKyBnZW9tX2hpc3RvZ3JhbShiaW53aWR0aCA9IC4wMSkgKyBmYWNldF93cmFwKCB+IHNfbnVtKQoKZmlsdF9sYWcgPC0gbGFnX2RhdGEgJT4lIGdyb3VwX2J5KHNfbnVtKSAlPiUgZmlsdGVyKGxvZ19sYWcgPCAobWVhbihsb2dfbGFnKSArIDEuNSpzZChsb2dfbGFnKSkgLCBsb2dfbGFnID4gKG1lYW4obG9nX2xhZykgLSAxLjUqc2QobG9nX2xhZykpKSAlPiUgdW5ncm91cCgpCgpnZ3Bsb3QoZmlsdF9sYWcsIGFlcyh4PWxhZykpICsgZ2VvbV9oaXN0b2dyYW0oYmlud2lkdGggPSAxMCkgICsgZmFjZXRfd3JhcCh+c19udW0pICsgeGxpbSg1MCwzNTApCgpnZ3Bsb3QoZmlsdF9sYWcsIGFlcyh4PWxvZ19sYWcpKSArIGdlb21faGlzdG9ncmFtKGJpbndpZHRoID0gLjAyKSArIGZhY2V0X3dyYXAofnNfbnVtKQoKYGBgCgpgYGB7cn0KIyMgd2l0aCBhIHJ1biwgcGxvdCB0aGUgbWVhbiBsYWcgCnJvdW5kX2ZpbHRfbGFnIDwtIGZpbHRfbGFnICU+JSAKICBncm91cF9ieShzX251bSwgbl90cmF2ZWxfc3RlcHMsIHRyYXZlbF9rZXlfY29uZCwgcm91bmQsIHBoYXNlKSAlPiUKICBzdW1tYXJpc2UodHJpYWxfbnVtID0gZmlyc3QodHJpYWxfbnVtKSwgbWVhbl9sYWcgPSBtZWRpYW4obGFnKSwgbWVhbl9sb2dfbGFnID0gbWVkaWFuKGxvZ19sYWcpKSAlPiUgbXV0YXRlKHJvdW5kX251bSA9IGFzLmludGVnZXIoYXMuY2hhcmFjdGVyKHJvdW5kKSkpCgpmb3IgKHMgaW4gMToyMil7CiAgCiAgIyBwbG90IG1lYW4gbGFnIGFzIGEgZnVuY3Rpb24gb2Ygcm91bmQgbnVtYmVyCiAgcCA8LSBnZ3Bsb3Qocm91bmRfZmlsdF9sYWcgJT4lIGZpbHRlcihzX251bSA9PSBzKSwgYWVzKHggPSByb3VuZF9udW0sIHkgPSBtZWFuX2xvZ19sYWcsIGNvbG9yID0gcGhhc2UpKSArIGdlb21fcG9pbnQoc2l6ZSA9IDIpICtmYWNldF9ncmlkKG5fdHJhdmVsX3N0ZXBzIH4gdHJhdmVsX2tleV9jb25kKSArIGdlb21fbGluZShzaXplID0gMS4yKSAgKyB4bGFiKCdoYXJ2ZXN0L3RyYXZlbCBudW1iZXInKSArIHlsYWIoJ21lYW4gbG9nIGxhZycpICsgZ2d0aXRsZShwYXN0ZSgnc3ViajogJywgcykpICsgdGhlbWVfbWluaW1hbCgpCiAgcGxvdChwKQogIAp9CiAKCmBgYApgYGB7cn0KbGlicmFyeShwbG90cml4KQpjb25kX3JvdW5kX2ZpbHRfbGFnIDwtIHJvdW5kX2ZpbHRfbGFnICU+JSAKICBncm91cF9ieShuX3RyYXZlbF9zdGVwcywgdHJhdmVsX2tleV9jb25kLCByb3VuZF9udW0sIHBoYXNlKSAlPiUKICBzdW1tYXJpc2UobWVhbl9sb2dfbGFnID0gbWVkaWFuKG1lYW5fbG9nX2xhZyksIHNkX2xhZyA9IHN0ZC5lcnJvcihtZWFuX2xvZ19sYWcpKQoKcCA8LSBnZ3Bsb3QoY29uZF9yb3VuZF9maWx0X2xhZywgYWVzKHggPSByb3VuZF9udW0sIHkgPSBtZWFuX2xvZ19sYWcsIGNvbG9yID0gcGhhc2UpKSArIGdlb21fcG9pbnQoc2l6ZSA9IDIpICtmYWNldF9ncmlkKG5fdHJhdmVsX3N0ZXBzIH4gdHJhdmVsX2tleV9jb25kKSArIGdlb21fbGluZShzaXplID0gMS4yKSArIGdlb21fbGluZXJhbmdlKGFlcyh5bWluID0gbWVhbl9sb2dfbGFnIC0gc2RfbGFnLCB5bWF4ID0gbWVhbl9sb2dfbGFnICsgc2RfbGFnKSkgKwogIHhsYWIoJ2hhcnZlc3QvdHJhdmVsIG51bWJlcicpICsgeWxhYignbWVhbiBsb2cgbGFnJykgKyBnZ3RpdGxlKCdncm91cDogJykgKyAgdGhlbWVfbWluaW1hbCgpCgpwbG90KHApCmBgYApgYGB7cn0KY29uZF9yb3VuZF9maWx0X2xhZwpgYGAKCgpgYGB7cn0KCnRyaWFsX2ZpbHRfbGFnIDwtIHJvdW5kX2ZpbHRfbGFnICU+JSAKICBncm91cF9ieShzX251bSwgbl90cmF2ZWxfc3RlcHMsIHRyYXZlbF9rZXlfY29uZCwgcGhhc2UpICU+JQogIHN1bW1hcmlzZSh0cmlhbF9udW0gPSBmaXJzdCh0cmlhbF9udW0pLCAKICAgICAgICAgICAgbWVhbl9sYWcgPSBtZWFuKG1lYW5fbGFnKSwgCiAgICAgICAgICAgIG1lYW5fbG9nX2xhZyA9IG1lYW4obWVhbl9sb2dfbGFnKSkgJT4lIHVuZ3JvdXAoKSAlPiUgbXV0YXRlKHNfbnVtID0gYXMuZmFjdG9yKHNfbnVtKSkKCmdncGxvdCh0cmlhbF9maWx0X2xhZywgCiAgICAgICBhZXMoeCA9IGZhY3RvcihuX3RyYXZlbF9zdGVwcyksIHkgPSBtZWFuX2xvZ19sYWcsIGZpbGwgPSBzX251bSkpICsgCiAgZ2VvbV9wb2ludChhZXMoY29sb3IgPSBzX251bSksIHNpemUgPSAyKSsgCiAgZ2VvbV9saW5lKGFlcyhncm91cCA9IHNfbnVtLCBjb2xvciA9IHNfbnVtKSwgc2l6ZSA9IDEuMikgKwogIGZhY2V0X2dyaWQocGhhc2UgfiB0cmF2ZWxfa2V5X2NvbmQpICsgdGhlbWVfbWluaW1hbCgpICsgeGxhYignIyB0cmF2ZWwgc3RlcHMnKSArIHlsYWIoJ2xvZyBsYWcgdGltZScpCgpnZ3Bsb3QodHJpYWxfZmlsdF9sYWcsIAogICAgICAgYWVzKHggPSB0cmF2ZWxfa2V5X2NvbmQsIHkgPSBtZWFuX2xhZywgZmlsbCA9IHNfbnVtKSkgKyAKICBnZW9tX3BvaW50KGFlcyhjb2xvciA9IHNfbnVtKSwgc2l6ZSA9IDIpKyAKICBnZW9tX2xpbmUoYWVzKGdyb3VwID0gc19udW0sIGNvbG9yID0gc19udW0pLCBzaXplID0gMS4yKSArCiAgZmFjZXRfZ3JpZChwaGFzZSB+IG5fdHJhdmVsX3N0ZXBzKSArIHRoZW1lX21pbmltYWwoKQpgYGAKCmBgYHtyfQoKCmdyb3VwX2xhZyA8LSB0cmlhbF9maWx0X2xhZyAlPiUgZ3JvdXBfYnkobl90cmF2ZWxfc3RlcHMsIHRyYXZlbF9rZXlfY29uZCwgcGhhc2UpICU+JSBzdW1tYXJpc2UoZ21fbGFnID0gbWVhbihtZWFuX2xhZyksIGdtbF9sYWcgPSBtZWFuKG1lYW5fbG9nX2xhZykpICU+JSB1bmdyb3VwKCkKCmdncGxvdChncm91cF9sYWcsIGFlcyh4ID0gZmFjdG9yKG5fdHJhdmVsX3N0ZXBzKSwgeSA9IGdtbF9sYWcsIGNvbG9yID0gdHJhdmVsX2tleV9jb25kKSkgKwogIGdlb21fbGluZShhZXMoZ3JvdXAgPSB0cmF2ZWxfa2V5X2NvbmQpLCBzaXplID0gMikgKyBmYWNldF93cmFwKH5waGFzZSkgKyBnZW9tX3BvaW50KHNpemUgPSAyKSAreGxhYignIyB0cmF2ZWwgc3RlcHMnKSArCiAgeWxhYignZ3JvdXAgbWVhbiBsb2cgbGFnJykgKyB0aGVtZV9taW5pbWFsKCkKCgoKYGBgCmBgYHtyfQpoYXJ2ZXN0X2xhZyA8LSByb3VuZF9maWx0X2xhZyAlPiUgdW5ncm91cCgpICU+JSBmaWx0ZXIocGhhc2UgPT0gIkhBUlZFU1QiKSAlPiUgbXV0YXRlKHN1YmogPSBhcy5mYWN0b3Ioc19udW0pKQpoZWFkKGhhcnZlc3RfbGFnKQoKaGxfbW9kZWwgPC0gbG1lcihzY2FsZShtZWFuX2xvZ19sYWcpIH4gIG5fdHJhdmVsX3N0ZXBzICsgdHJhdmVsX2tleV9jb25kICsgcm91bmRfbnVtICArICgxICsgbl90cmF2ZWxfc3RlcHMrIHRyYXZlbF9rZXlfY29uZCArIHJvdW5kX251bSB8c3ViaiksIGRhdGEgPSBoYXJ2ZXN0X2xhZywgY29udHJvbCA9IGxtZXJDb250cm9sKG9wdGltaXplciA9Im9wdGlteCIsb3B0Q3RybD1saXN0KG1ldGhvZD0nbmxtaW5iJykpKQoKc3VtbWFyeShobF9tb2RlbCkKYGBgCgoKCmBgYHtyfQpsaWJyYXJ5KGJyb29tLm1peGVkKQp0ZCA8LSB0aWR5KGhsX21vZGVsLCBjb25mLmludCA9IFRSVUUpCnRkc29tZSA8LSB0ZCAlPiUgZmlsdGVyKHRlcm0gPT0gJ25fdHJhdmVsX3N0ZXBzJyB8IHRlcm0gPT0gJ3RyYXZlbF9rZXlfY29uZEhBUkQnIHwgdGVybSA9PSAncm91bmRfbnVtJykKCnRkc29tZSA8LSB0ZHNvbWUgJT4lIG11dGF0ZShwYXJhbWV0ZXIgPSBjYXNlX3doZW4odGVybSA9PSAnbl90cmF2ZWxfc3RlcHMnIH4gJ3RyYXZlbCBzdGVwJywgdGVybSA9PSAndHJhdmVsX2tleV9jb25kSEFSRCcgfiAndHJhdmVsIGVmZm9ydCBjb25kJywgdGVybSA9PSAncm91bmRfbnVtJyB+ICdyb3VuZCBudW1iZXInKSkKCmdncGxvdCh0ZHNvbWUsIGFlcyhlc3RpbWF0ZSxwYXJhbWV0ZXIsY29sb3IgPSBwYXJhbWV0ZXIpKSArIAogIGdlb21fcG9pbnQoc2l6ZSA9IDIpICArCmdlb21fZXJyb3JiYXJoKGFlcyh4bWluID0gY29uZi5sb3csIHhtYXggPSBjb25mLmhpZ2gpLCBzaXplID0gMS4yKSArIHRoZW1lKGxlZ2VuZC5wb3NpdGlvbiA9ICJub25lIikgKyB0aGVtZV9taW5pbWFsKCkgKyB5bGFiKCdmaXhlZCBlZmZlY3QnKSArIHhsYWIoJ2VzdGltYXRlICgjIHBvaW50cyknKSArIGdndGl0bGUoJ2xpbmVhciBtaXhlZCBlZmZlY3RzIG1vZGVsIC0gZWZmZWN0IG9uIGxhZycpICsgdGhlbWUoYXhpcy50ZXh0LnkgPSBlbGVtZW50X3RleHQoYW5nbGUgPSA5MCwgaGp1c3QgPSAuNSkpCgpgYGAKCg==